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RAY TRACING FOR IONOSPHERIC AND TROPOSPHERIC PROPAGATION STUDIES 



SUMMARY 

Adapted for a digital computer, conventional ray tracing becomes a 
potentially useful technique for studying several aspects of long-range wave- 
propagation through stratified media. The development of a ray-tracing pro- 
gramme and its application to both medium-frequency sky-wave propagation 
through the ionosphere and ultra-high-frequency (u.h.f.) propagation in the 
troposphere are described in this report. 



1. INTRODUCTION 

In propagation through inhomogenous or strati- 
fied media, methods based on a rigorous wave- 
theory solution are usually cumbersome and their 
practical application lengthy, even if it can be 
assumed that the model representing the medium is 
accurate and realistic enough to justify such an 
elaborate treatment. On the other hand, a geomet- 
rical-optics approach, although only an approxi- 
mation, is relatively simple in conception and 
application. In fact, it has been shown in many" 
instances to lead to a surprisingly good estimate of 
the general behaviour of propagated waves. Further- 
more, ray tracing based on geometrical -optics 
principles is readily adapted to evaluation by com- 
puter, because of its simple iterative nature. 

It is not intended to discuss here the degree of 
validity of a ray-theory approach to long-range 
propagation. We assume from the outset that the 
variations in the refractive index of the medium 
with height are sufficiently slow compared with the 
wavelength and, in the application to long-range 
tropospheric propagation that the atmospheric duct 
guiding the radiation is sufficiently extensive, 
horizontally and vertically. Even so, the ray- 
tracing interpretation of the well-known ducting 
phenomenon in the troposphere (at u.h.f. for example) 
can produce a complicated pattern of rays indicating 
the possibility of strong fields at certain points in 
the propagation range. (An example of this can be 
seen in Fig. 5.) This interpretation cannot readily 
give quantitative results when strong ducting is 
present but can indicate whether a certain refractive 
index profile would be expected to favour enhanced 
fields at great distances. 

For application to ionospheric propagation, the 
method developed is expected to be valid for both 
the medium-frequency and high-frequency broadcast 
bands. At frequencies less than the gyromagnetic 
frequency (approximately 1 MHz) the method should 



not be used for waves launched with angles of 
elevation greater than about 30°. This restriction 
can be relaxed for high-frequency propagation. 

One application of ray tracing at medium fre- 
quencies is the determination of absorption, and it 
may be important in some applications (e.g. in the 
calculation of ionospheric cross-modulation) to 
know not only the total absorption for the path but 
also the relative amounts of absorption that occur 
at different points along the path. The method of 
calculation is therefore extended to include absorp- 
tion, as discussed in Section 4.3. 

The method is based on the well established 
principle that the wave-normal obeys Snell's law of 
refraction and on the assumption that the propaga- 
tion characteristics of the medium are spherically 
stratified, concentric with a spherical Earth. Snell's 
law, applied to spherically curved layers takes the 
form (Fig. 1). 
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Fig. 1 - Snell's law applied to a spherically 
stratified medium 



fiTsin (fi = constant 



(1) 



where, at a given point in the medium, P in Fig. 1 
say, 

fi = the real part of the refractive index of the 
medium at P 



r = the radius of curvature 
index sphere through P 



of the constant- 



and (ji = the angle between the wave-normal and the 
normal to the constant-index sphere at P, 
i.e. the angle of incidence at P. 
Clearly, because of the initial assumptions, r will 
be a linear function of the vertical height h of the 
ray above the ground and the refractive index will 
be- a continuous function of h. Thus, by virtue of 
Equation (1), the direction of the ray will be a 
function of h, and if this is computed for a sufficient 
number of points the ray trajectory is established. 

The application to tropospheric propagation at 
u.h.f. is simplified by the fact that the refractive 
index is everywhere real and independent of the 
Earth's magnetic field, since the numbers of free 
electrons and gas ions are insignificant. Here, \i is 
a function of h only. 

Waves penetrating into the magneto-ionic 
plasma associated with the ionosphere, however, 
are influenced by the Earth's magnetic field in a 
complicated manner^-^. The incident wave is split 
into two components, the ordinary and extraordinary 
waves, which travel generally in different direc- 
tions, and energy is absorbed from each of them by 
the medium. Furthermore, at any point the wave 
energy in each component is propagated in a slightly 
different direction to the direction of the wave- 
normal and is not confined to the initial plane of 
incidence. Although the latter effect can be allowed 
for in the ray-tracing calculations, it is neglected 
in the treatment described below because, for the 
applications of the method envisaged, it is believed 
that the resulting error will be small. 

The complex refractive index of the magneto- 
ionic medium for the ordinary and extraordinary 
components is taken to be that predicted by the 
Appleton-Hartree formulae.* These formulae show 
that the real part of the complex refractive index is 
a function of several variables, one of which is the 
angle between the ray direction and that of the 
Earth's magnetic field. This fact complicates the 
situation somewhat because Snell's law (Equation 
(1)) cannot then be solved explicitly. In the step- 
by-step procedure outlined below, however, this 
difficulty can be largely overcome. 



These formulae are fully discussed in Ref. 1. 
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Fig. 2 - The ray trajectory: co-ordinates and 
notation 

2. THE RAY TRAJECTORY 

Fig. 2 shows the trajectory of a typical ray 
launched from a point near the Earth's surface with 
an angle of elevation a (take-off angle). Using 
polar co-ordinates (r, ?/) with the origin at the centre 
of the Earth and the diameter through the transmitter 
(77 = 0) for angular reference, the trajectory (which 
lies in a great circle plane) is described by r= F(77). 
Let Tg be the radius of the Earth and h the vertical 
height of a point on the ray above the ground, so 
that r =r^+h = F(7/). 

Consecutive points on the trajectory which are 
separated by a small angular interval, e say, are 
expected to be well-related by the first few terms of 
an infinite Taylor series. The series relation is: 



dF f^d=F e^^d^'F 

F(77 + €) = F(77) +«-_ + _—-+ - + 

ari 
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(2) 



Since the trajectories pertaining to wave-propagation 
problems are usually reasonably smooth and well- 
behaved, it is estimated that for a ground range 
interval A(= r^e see Fig. 2) of a few kilometres the 
first three terms of the Taylor series will suffice. 

Let r^^ 1 be the height the ray r = F(t]) has 
reached after the (k + 1)**^ interval; then from 
Equation (2) 



'■/c + 1 = ^fc + e 



\dr,j^^ 2\^rfj, 



(3) 



If /xo, To and 4>o are the initial values of these 
variables at the beginning of the trajectory we have 
from Equation (1), \jis'\n oS = ftoTosin <i^= constant. 



Making the substitution sin 00 - cosa and rearrang- 
ing gives 



cosec <ji 



°rjv7 



seca 



(4) 



where (Fig. 2) 

a= the tal<e-off angie 

To = the value of r for the transmitting point 

^Lo = the refractive index at r = tq. 
Now it may be shown that 
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(5) 



and, by differentiating Equation (1) to obtain 
(d^/dr), that 
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cosec^(^;(^ - 1 > (6) 



of radiation by an abnormal combination of refrac- 
tive-index gradients in the lower regions of the 
troposphere. The conditions which are conducive 
to such propagation phenomena are readily explored 
by ray-tracing techniques. It must be pointed out, 
however, that if the atmospheric duct is regarded as 
a waveguide, the possibility exists of the latter 
being beyond cut-off for the wavelength considered 
If the effective height of the guide is too small. 
This possibility is not indicated by a ray-tracing 
approach which must, therefore, be used with some 
caution in this application. 

The refractive index gradients of the air layers 
are very weak so that propagation much beyond the 
normal horizon occurs only for waves launched with 
very small take-off angles, say -1°<a<1° and in 
Equation (7) sec^a can be replaced by 1 + cx^. For 
similar reasons the maximum height of any signifi- 
cant ray is not likely to exceed 2 km above the 
ground. The refractive index, /x, of the air in the 
troposphere is always near unity and, therefore, it 
is more convenient to specify the index in i\'-units, 
i.e., as N = (n - 1)10^. 

Bearing in mind these limitations on the range 
of the variables, Equation (7) can be reduced to 
read 



/rfc, 1 =^k +A0 -08463% 0-314 



0-0846a%0-314(/ifc - hj] + 2^^ -Nj) 
1 +6-371 J 

L V'^^ A. 



whereupon, from Equations (4), (5) and (6), Equation 
(3) becomes 



'■fc+ 1 ° Tfe + rfcf 






+ 0-0785A 
where 



(8) 



hjf^ = height of ray above the ground in metres 
at the ground range feA, fe = 1, 2, 3, etc. 

/jj = height of transmitting aerial in metres 

Njt- = refractive indexat the height hj^ in N-units 

/V J = refractive index at the height hj jn I'V- 
unlts 



{^J 



a = take-off angle in minutes of arc 



sec a - 1 , 






Equation (7) is the basic iterative relation which is 
used to compute the ray trajectory. Below, further 
reductions of Equation (7) are developed which are 
better adapted to the characteristics of the variables 
in the particular application. 



3. LONG-RANGE U.H.F. PROPAGATION IN THE 
TROPOSPHERE 

3.1. General Formulation 

It is well known that one mechanism for long- 
range propagation at u.h.f. is the ducting or guiding 



gradient of the refractive- index at the 
height h^. (N-units/ metre). 

A = selected ground-range interval in kilo- 
metres. 

The Earth's radius* is taken to be 6-37 x 10'' metres. 

3.2. Sign Determination 

Measurements of the refractive index of air 
often reveal a very irregular distribution with height 
and the gradient d.N/d/i may be positive or negative, 

* Note that a smooth, spherical Eorth was assumed at 
the outset. 



although usually the latter. It will be seen from 
Equation (8) that the sign of the A'^ term is uniquely 
determined by the magnitude and sign of the gradient 
dN/6h and, therefore, by the input data. The 
second term in Equation (8), however, contains a 
square-root factor and whether the positive or 
negative root is taken depends on whether the ray 
height is increasing or decreasing respectively. 
The ray-tracing programme must contain instructions 
to determine the correct sign for this term at each 
step in the calculation. 

The following procedure is found to be satisfactory: 

(a) Initial step: choose the positive sign if a is 
positive and the negative sign if a is negative. 
If a = choose the same sign as the A^ term. 

(b) Subsequent steps: Fom hj^+t - h using the 
sign determined by (a), and change the sign 
whenever |/tfc+ i - /z/c |^|A^ term|. (This rule 
determines when the ray height has passed 
through a maximum or a minimum.) 

3.3. Input Data 

The refractive- index profile is presented as a 
sequence of height-index pairs of values (height in 
metres, index in N-units). The height values are 
not in general equispaced but are given in increas- 
ing order. The index N^ corresponding to the height 
hj^ is obtained by first-order interpolation of the 
given pair-values. The gradient (6N/6h\ is 
assumed to be constant between successive pairs 
of values: thus if /j;^ lies between h^ and /ij + i the 
gradient is taken to be (Nj + i - Nj)/(/ij+ i - ^j). 

The range interval A, the transmitting-aerial 
height, and the take-off angles required constitute 
the rest of the input data. 

3.4. Output Data 

The data printed out depends largely on the 
form of the ray-trace investigation. A useful ex- 
ploratory arrangement is to compute the trajectory 
for each of a fan of rays with take-off angles centred 
on a = 0, e.g. -0°10' to +0°10' in 0-5' steps. The 
calculation is stopped when the ray either hits the 
earth or 'escapes', i.e. when the ray height exceeds 
the maximum height of the profile input data. The 
ray height at each step (or at intervals of several 
steps) is printed out if a detailed study of the 
trajectory is required. On other occasions, printing 
out only those portions of the trajectory near the 
turning points will suffice. 



(fx-lx) =i-x;i-iz 



Crude ray-focusing occurs with some index 
profiles and an estimate of the convergence gain in 
these caustic regions can be obtained by printing 
out the ray heights, for each of the fan of rays, at 
a specific ground range. The number of rays per 
unit height relative to that expected under free- 
space conditions (straight-line rays) gives ameasure 
of the power gain. 
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Fig. 3 - Oblique ionospheric propagation 



4. M.F. AND H.F. SKY-WAVE PROPAGATION VIA 
THE IONOSPHERE 

In this problem straight-line propagation through 
the tropospheric region is assumed but refraction 
occurs on entering the ionospheric region. Waves 
entering at sufficiently oblique angles of incidence 
(Fig. 3) are deviated back to Earth in the E layer 
after passing through the D layer, giving rise to 
long-range propagation. Further, wave energy is 
lost to the magneto-ionic medium and this absorp- 
tion loss may differ considerably in magnitude for 
the ordinary and extraordinary waves respectively. 
Both the ray trajectory and the absorption loss can 
be determined, for each component, if the variation 
of the complex refractive index with height is 
known. Each component is treated separately and 
it is assumed that an ordinary wave remains ordinary 
while traversing the ionosphere. This is justified 
in the absence of large changes of refractive index 
within height intervals comparable with a wave- 
length. If sudden changes are present, propagation 
will involve partial reflections and simple ray- 
tracing techniques will not be accurate. 

The calculations do not include polarization 
coupling losses, i.e. do not allow for the fraction 
of the power in the wave radiated from the trans- 
mitting aerial (usually linearly polarized) that 
excites the ordinary or extraordinary mode on enter- 
ing the ionosphere. 

4.1. Complex Refractive Index 

The Appleton-Hartree formula for the complex 
refractive index (ii-'jx) is given by 



y^sin^^ 
■2(1-X-jZ) 



y*sin*(9 
4(1-X-jZ)^ 



+ Y'cos'f? 



'! 



(9) 



where (i = real part of the refractive index 

X= imaginary part of the refractive Index. 
(N.B. y = If Z =0 or X = 0). 

X = 81 E/f ^; E is the free-electron density 
in number/cm° and f is the wave fre- 
quency in kHz. 

Y = Gyromagnetic frequency divided by the 
wave frequency. 

= angle between geomagnetic field and 
ray direction. 

2?7Z = collision frequency divided by the 
wave frequency. 

(+) =+ sign for ordinary wave, - sign for 
extraordinary wave. 

The angle (9 in Equation (9) Is calculated from the 
trigonometric formula 

cos = cos vV sine/) cosy - sin fA cos ci (10) 

where i/f = angle of magnetic dip: range -90° to 
+ 90° (positive in northern latitudes). 
N.B. the dip angle may vary significantly 
even over that part of the trajectory 
which is through the ionosphere. Usual- 
ly, however, it is sufficient to take its 
value at the mid-point of the path. 

4> = angle between the upwards vertical and 
the ray direction: range 0° to 180°. 

y = ray magnetic bearing: range 0° to 360° 
east of magnetic north. N.B. In addition 
to variations due to changes of magnetic 
declination along the path the ray bear- 
ing will be subject to a systematic vari- 
ation but, except possibly for very high 
latitudes, variations from the initial 
value may be neglected in the present 
application. 

The variations of electron density and collision 
frequency with height form the main input data, to- 
gether with the other necessary physical constants. 
It should be noted that the angle t^, which is re- 
quired to evaluate Equation (10), is a function of 
the ray direction (see Equation (1)) and cannot 
therefore be obtained directly from the input data. 

4.2. Calculation Procedure 

As in the previous (tropospheric) application, 
a ground range interval is selected and the ray 
height calculated at the end of each interval using 
the basic iterative relation of Equation (7). It is 
convenient to re-state the latter in the form 



{' 



P;t4i =Pfe +P;cf (P/t f"c sec a-1 +Pfc-" 



2 + 



where p^ = relative height (jy^lr^ = 1 + h^r ^ 
(re = 6-37 X 10® metres) 

i = angular range interval (range correspond- 
ing to height p^ is kr^i). 

n^ = real part of refractive index at the rela- 
tive height p^. 

a = take-off angle (transmitting aerial assum- 
ed at ground level). 




Fig. -i 



(b) 
Calculation procedure 



(a) points at entry into ionosphere 
(b) points at the turn-over 

Consider (Fig. 4(a)) points on that portion of the 
trajectory where it enters the ionosphere, i.e. at 
the height where X in Equation (9) first becomes 
significantly different from zero. Up to this height 
the refractive index is assumed real and equal to 
unity so that no absorption or refraction has occur- 
red. In Fig. 4(a), p^ marks the first step into the 
ionosphere and ^^ and y^. are calculated using 
Equations (9) and (10) with the approximation 



qS;^ « arc tan i 



Pk-Pk- 



(12) 



Next, the index gradient (dfi/dp)fc is obtained from 
the approximation 




Pk-Pk-i 



(13) 



Finally, Pk+i '^ calculated from Equation (11) 
the procedure repeated for subsequent steps.* 



and 



* The sign of the square root factor in Equation 11 is 
determined by the rules given in Section 3.2. 
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It will be seen that the approximation in Equa- 
tion (13) can become rather indeterminate at the 
turn-over point of the trajectory (Fig. 4{b)) since 
both numerator and denominator may then approach 
zero. Furthermore, at this point, <;6fc passes through 
7t/2 and the second term on the right-hand side of 
Equation (11) vanishes, so that the third term con- 
taining the index gradient (dff/dp)k is here very 
significant. In order to avoid a possible distortion 
of the ray trajectory at the turnover the index 
gradient is calculated from the partial derivative 
of the Appleton-Hartree formula, w/ith respect to 
the relative height p, w/henever \pic-Pk.i\<g' where 
g is an arbitrarily chosen small quantity. The co- 
efficient (d/i/dp)fc in Equation (11) is then replaced 
by (dji/dp)f^, the formula for which is given in the 
Appendix. 

4.3. Absorption 

An important part of the application to m.f. 
propagation is the ability to assess the absorption 
losses and the way these vary with height, range 
and frequency. As previously mentioned the imag- 
inary part x> of the complex refractive index, is 
calculated at the end of each interval simultaneous- 
ly with the real part fi. 

Let Xjt. 1 and Xfc be successively computed 
values of the imaginary part. The average value is 
(Xk.i + X)^)f2 and it may be shown that the average 
rate of energy loss over this interval is 

0-09083 { (Xfc. 1+ Vfc) decibels/km 

(f = wave frequency in kHz) 

The length of the interval along the ray trajectory 
is approximately 



hk-Pk-A 



2 2 

+ Pk. 



Vi 



so that the absorption loss, L^ ^, ^ say, is given 
by 



medium-frequency propagation.* A programme sub- 
routine carries out a first-order interpolation of 
these data to obtain the values at intermediate 
heights. A ray is deemed to have 'escaped' and 
further calculation is stopped if the calculated ray 
height exceeds the highest value entered for the 
profile. 

Other essential input data comprise the take- 
off angle, ray bearing, magnetic-dip angle, gyro- 
frequency, wave frequency and range interval. 
Separate calculations are necessary for the ordinary 
and extraordinary rays and the required case is 
selected by entering the reference number +1 or -1 
respectively. 

4.5. Output Data 

The range, ray height and cumulative absorp- 
tion loss are printed-out for each step in the ionos- 
phere. If the ray returns to earth the total range is 
calculated and printed-out, otherwise 'ray escapes' 
is printed. 

4.6. Additional Programme Features 

One application of ray tracing is to estimate 
the ionospheric absorption loss for propagation over 
a specific link. In this case the geographical 
locations of the transmitter and receiver are known 
and interest centres on ray traces within a narrow 
range of take-off angles. A sub-routine has been 
included in the main programme to deal with this 
requirement. 

The latitude and longitude of both terminal 
points are entered as input data from which the 
total range and magnetic bearing required are 
calculated, using straightforward trigonometric 
formulae. Further, again using standard formulae, 
an estimate is made of the nominal take-off angle 
required (a typical virtual height of about 100 km 
is assumed for E-layer reflection). Three ray traces 
are then computed for the take-off angles (aj-1)°, 
cxj° (aj+1)° where a- is the nominal angle required. 



Lfc,fc. 1= 0.09083 /•r,(X^.^ + Xfc) 
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In the ray-tracing programme the running-total 
absorption loss 

1 
is calculated and printed out together with the 
range and ray height at the end of each interval. 

4.4. Input Data 

The electron-density and collision-frequency 
profiles of the mode! ionosphere are entered as a 
sequence of values at intervals of one kilometre, 
normally extending from about 60km to 110km for 



The corresponding ground ranges so obtained are 
generally found to straddle the exact range required 
and, by subsequent graphical interpolation or 
inspection, it is easy to see how the absorption 
loss behaves in the region of interest. Also, from 
the increments of ground range between the succes- 
sive rays, the actual convergence gain can be de- 
duced. However, it has been found that this differs 
little from the convergence gain produced by a 
spherical mirror at the virtual height. 

For high-frequency applications, where the reflection 
heights ore much greater, the height intervals at which 
data are entered and the upper limit con be increased 
accordingly. 
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Fig. 5 - Ray set in the troposphere: illustrating ducting phenomena 
a = take-off angle in minutes of arc 
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5. ILLUSTRATIVE RESULTS 

Fig. 5 shows a ray set obtained in the tropos- 
pheric application when the refractive-index grad- 
ients of the air layers close to the ground are con- 
ducive to the formation of elevated ducts. The 
mechanism of long-range propagation from an elev- 
ated transmitting aerial under these conditions is 
clearly demonstrated. The number labelling each 
ray in the figure is the take-off angle in minutes 
of arc. The computing time (including print-out) 
using the Elliot 803B is approximately one minute 
per ray. 

Fig. 6 shows the trajectories of several ord- 
inary rays through the ionosphere with a fixed take- 
off angle but varying wave frequencies. An elec- 
tron-density profile for the ionosphere was assumed 
which corresponded to average night-time conditions 
found in Europe. The greater penetration and range 
of the higher-frequency waves is evident. The 
asymmetry of the trajectories is due to the difference 
in the angle between the ray direction and Earth's 
magnetic field for the upward and downward paths. 

Fig. 7 shows a plot of ray height versus cum- 
ulative absorption loss for two ordinary rays as 
they progress through the ionosphere. The rays 
have the same wave frequencies and are launched 
by a transmitter in an equatorial location with the 
same take-off angle, and the ground ranges are 
approximately equal. The solid-line upper curve 
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Fig. 6 - Ordinary rays in the ionosphere: 
effect of wave frequency 

refers to east-west propagation while the dashed- 
line curve refers to north-south propagation. The 
lower absorption loss for north-south propagation in 
in equatorial regions is clearly indicated. 



6. CONCLUSIONS 



fgyro =0 9 MHz 
f =1 5 MHz 

(2 = 9° 

ground range = 1100 km 
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Fig. 7 - Ionospheric absorption losses (ordinary 
waves): variation with height 

East-west propagation at the equator 



A simple computer programme for conventional 
ray tracing, such as that outlined in this report, 
is certainly a useful tool for exploratory research 
into several aspects of propagation phenomena. It 
is attractive because of the rapidity with which 
results can be obtained, the control of the para- 
meters which the investigator can exercise, and the 
potential adaptability to media with a wide range 
of characteristics. Extension of the present pro- 
gramme to include other more complicated fomis of 
layer stratification should not prove difficult. 

It must be emphasized that, even for those 
conditions where a geometrical ray-theory is a 
valid approximation to a full wave treatment, the 
realism of the results obtained is governed largely 
by the accuracy with which the physical character- 
istics of the stratified media can be ascertained by 
measurement. 
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APPENDIX 



Calculation of Refractive-Index Gradient 



Normally, the index gradient is most easily 
calculated from the successively computed values 
of n. Thus the differential coefficient (d^/dp)jj 
required in the trajectory equation Is obtained by 
calculating (p^k-i^k- i}^(Pk~Pk- ■)• ^^ *^6 top o) the 
trajectory, however, it may happen that p^ becomes 
very close or even equal to pj.. i and the gradient 
calculated in this manner is likely to be in error, 
leading to distortion of the ray trajectory and errors 
In the absorption loss and ground range. In this 
event, i.e. when Pk^Pk-j' it is better to calculate 
the gradient from the derivative of the Appleton- 
Hartree formula, with respect to the relative height 
p. As a first approximation one may use the partial 
derivative (dfi/dp)i^. The approximation will be good 
if the rate of change of the ray direction (which, 
unfortunately, is a maximum at the turn-over point) 
does not significantly affect the index gradient. 

In order to simplify the calculation, it is assum- 
ed that the collision frequency does not affect the 
ray path near the turn-over point so that Z can be 
set equal to zero in the Appleton-Hartree formula. 
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Thus from 




2{1-X) 

By definition X = 8lE/f^, therefore 

dX X dE 

dp E dp 

Hence substituting and inserting subscripts the 
partial derivative is: 



(l-f'fc) (fifc(1-f^/?) 



2(^fe Ek 



^-Xu 



Xfe-(i-/^fc) 



X^+{U^-^){^-n„l 
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